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The description of general relativistic radiation hydrodynamics in spherical symmetry is pre- 
sented in natural coordinate choices. For hydrodynamics, comoving coordinates are chosen, and the 
momentum phase space for the radiation particles is described in comoving frame four-momenta. We 
also investigate a description of the momentum phase space in terms of the particle impact param- 
eter and energy at infinity and derive a simple approximation to the general relativistic Boltzmann 
equation. Further developed are, however, the exact equations in comoving coordinates, because 
the description of the interaction between matter and radiation particles is best described in the 
closely related orthonormal basis comoving with the fluid elements. We achieve a conservative and 
, concise formulation of radiation hydrodynamics that is well suited for numerical implementation by 

a variety of methods. The contribution of radiation to the general relativistic jump conditions at 
shock fronts is discussed, and artificial viscosity is consistently included in the derivations in order 
to support approaches relying on this option. 
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I. INTRODUCTION 



J> ■ Astro-physical knowledge is gained by the observation of luminous objects and the search for scenarios that explain 
the observations based on accepted physics. Often, numerical simulations are required to determine if a detailed 
scenario really produces the given observation. In supernova theory, for example, the transition from analytical con- 
siderations of core collapse and explosion (Colgate and Johnson M) to numerical simulations (Colgate and White 
was made very early, when computers became available. General relativistic simulations of core collapse followed 
immediately (May and White J| ) . These simulations were based on the newly derived Einstein equations in spherical 
symmetry (Misner and Sharp |^f) that served as a basis for many later investigations, including also the present paper. 
In comoving coordinates, general relativistic radiation transport was first formulated by Lindquist |^]. Here, we refor- 
mulate general relativistic radiation hydrodynamics in light of modern numerical algorithms requiring conservation 
laws, and hope to reveal some of the beauty in the spherically symmetric case that may have remained hidden. 
However, comoving coordinates are only one possible choice out of a variety of 3+1 decompositions enabled by the 
+3 . covariance of general relativity in four-dimensional space-time (Arnowitt, Deser, and Misner ||; Smarr and York H). 
The Boltzmann transport equation is usually split into a left-hand side and a right-hand side. The left-hand side is 
the directional derivative of the distribution function along trajectories of free particle propagation. This derivative is 



equated to the changes in the distribution function due to collisions. Thus, the right-hand side accounts for changes 
in the radiation particle distribution function owing to particles that are created, annihilated, or scattered into new 
; ^ ■ states in the momentum phase space. 

After the choices of a 3+1 decomposition and a basis in the momentum phase space have been made, the directional 
derivative along the phase flow can be expressed in terms of partial derivatives of the distribution function with respect 
to the coordinates. The complexity of the left-hand side as well as the complexity of the collision term depends on 
the coordinate choice. A general discussion of radiation transport in spherically symmetric 3+1 decomposition has 
been provided by Mezzacappa and Matzner ||. Maximal slicing was chosen and the particle distribution function 
was described by the four-momenta measured in the frame of an observer comoving with the matter. This choice 
eases the evaluation of angle dependent cross sections in the collision term (Mihalas Mezzacappa and Bruenn 
fiof). Schinder and Bludman Jll[] chose polar slicing in the 3+1 decomposition and a tangent ray approach in the 
description of the momentum phase space. This choice avoids partial derivatives with respect to the momentum space 
variables on the left-hand side of the Boltzmann equation. 

In this work we explore the most natural 3+1 decomposition, one that is not enforced by an artificial external slicing 
condition. The time slices in the orthogonal comoving coordinates used by Misner and Sharp [Q are attached to the 
dynamical motion of the matter. As for the basis of the momentum space we first investigate the description in the 
comoving frame. In another part, we also develop a formulation in the spirit of the tangent ray approach and draw 
the connections between these two choices. But first of all, we continue with a motivation of the comoving coordinate 
choice. 
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The most obvious advantage is already contained in the word comoving. A given coordinate interval moves with 
the matter and adjusts to the location where the matter is. Comoving coordinates ease the description of the 
internal physical state of fluid elements because numerically difficult advection terms do not enter the hydrodynamics 
equations as in other coordinate systems. Gravitational collapse in spherical symmetry is a natural application of 
comoving coordinates. With regard to radiation transport, the comoving frame, which is preferred for the evaluation 
of the integrals over the angle dependent emission, absorption, and scattering kernels, is collinear with the comoving 
coordinates such that Lorentz transformations are avoided in the collision term. 

On the other hand, the orthogonal comoving coordinates show an important drawback in general relativity. After 
the formation of a black hole, observers may fall within finite time into the physical singularity at the symmetry 
center. In addition, a coordinate singularity forms at the Schwarzschild horizon. Simulations of the outer regions 
cannot proceed because of these singularities in the computational domain. 

To circumvent this difficulty, interest focused on singularity avoiding time slicing. For example, simulations of 
Wilson jl2j and Shapiro and Teukolsky (T^| were carried out in maximal slicing. The collapse to a black hole can 
be followed beyond the appearence of trapped surfaces in these coordinates. At late times however, the grid is 
"sucked down" the black hole, leading to unsatisfactory resolution in the physically interesting region outside the 
event horizon. Polar slicing shows even stronger singularity avoidance (Bardeen and Piran PH). The computational 
domain stays outside the apparent horizon for all times. This time slicing was implemented by Romero et al. |l5| l 
in orthogonal coordinates. Schinder et al. |ll| found an interesting compromise by introducing nonzero shift vectors 
to get comoving coordinates in polar slicing. Orthogonality of the coordinates, however, had to be abandoned. The 
idea of using observer time coordinates suggested by Hernandez and Misner |lj]] has been used by Miller and Motta 
jl8| and Baumgarte et al. in a singularity avoiding code. Although coordinate and physical singularities are not 
encountered in singularity avoiding schemes, the coefficients of the metric can increase without limits for late times 
(Petrich et al. [^0) ) . Alternatively, the computational domain can also be limited to the region outside of the event 
horizon by the choice of appropriate shift vectors (Liebendorfer and Thielemann pl| ) . A free choice of time slicing is 
then reestablished even in the presence of singularities. 

One argument against the use of comoving orthogonal coordinates is the difficulty to extend them to multiple spatial 
dimensions because of grid entangling. However, one-dimensional simulations offer the opportunity to include general 
relativity, exact Boltzmann radiation transport, and sophisticated physics for emission, absorption, and scattering at 
once. While the individual pieces are well represented in the literature, dynamical simulations with all ingredients are 
very difficult and about to become current state of the art. 

Many spherically symmetric simulations of compact objects were based on the comoving orthogonal coordinates 
of Misner and Sharp Q| . Finite difference schemes were constructed by May and White j|] , Van Riper |^2| , Bruenn 



J3|, Rezzolla and Miller 24 1, and Swesty p5[ . An approximate Riemann solver was constructed by Yamada [ p6[ 
However, the appearance of the time dependent metric in the general relativistic equations prevented an immediate 
application of numerical methods developed for conservatively formulated hydrodynamics. A formulation of the 
dynamics in terms of conservation laws has several benefits (i) Fundamental conservation laws are also valid at 
discontinuities and allow an accurate numerical solution - e.g., for the propagation of shock waves. This is not the 
case with arbitrary finite difference approximations, (ii) The integration of a conserved quantity over adjacent fluid 
elements does not depend on the fluxes at the enclosed boundaries. Thus, discretization errors have mainly local 
influence. This advantage is important for the implicit solution of problems involving scale differences of many orders 
of magnitude. (Hi) The integration over the whole computational domain of the single-fluid-element conservation laws 
leads to the total conservation of fundamental physical quantities, independent of the resolution in the conservative 
finite difference representation, (iv) An overwhelming variety of generalized investigations and numerical methods for 
the solution of hyperbolic conservation laws already exists (see, e.g., Davis [EtJ and references therein). In analogy to 



the work of Romero et al. [|15| in polar slicing, Liebendorfer and Thielemann |21| formulated conservative equations 
for hydrodynamics in the comoving frame. This approach is extended here to include radiation transport. The 
description (with the omission of neutrino back reaction to the fluid metric) has successfully been used in a general 
relativistic simulation of stellar core collapse, bounce and postbounce evolution, based on multigroup Boltzmann 
neutrino transport (Liebendorfer et al. papffl). 



II. CONSERVATIVE EINSTEIN EQUATIONS 

The left-hand side of the Einstein field equation, the Einstein tensor, is determined by the metric. Following Misner 
and Sharp B , we select in spherical symmetry 
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ds 2 = -a 2 dt 2 + - da 2 +r 2 (d$ 2 + sin 2 tidtp 2 ) , (1) 
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where r is the areal radius and a is a label corresponding to an enclosed rest mass (the prime denotes a derivative 
with respect to a: r' — dr/da). The proper time lapse of an observer attached to the motion of rest mass is related 
to the coordinate time dt by the lapse function a. The angles i? and ip describe a two-sphere. 

The right-hand side of the Einstein equations is given by the fluid- and radiation stress-energy tensor, T, which, in 
a comoving orthonormal basis, has the components 



T" = p(l + e) 



rpaa rjvd'd r j 1{ P l P p 



(2) 



The total energy is expressed in terms of the rest mass density, p, and the specific energy, e. The isotropic pressure 
is denoted by p, and radial net energy transport is included by the nondiagonal component q. Note our convention 
that energy density and pressure contain the contribution from the radiation field as well as the matter. 

In the Appendix, we rederive a system of equations equivalent to the Einstein equations in spherical symmetry, 
closely following the guideline of Misner et al. ||(J. The only difference to previous work is the omission of any 
approximations in the derivation. Moreover, the Appendix gives many useful relations for what follows. The final 
system of exact equations reads 
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The velocity u is equivalent to the r component of the fluid four-velocity as observed from a frame at constant 
areal radius (May and White ||). In the special relativistic limit, T = y/l + u 2 — 2m/r becomes the Lorentz factor 
corresponding to the boost between the inertial and the comoving observers. The gravitational mass m is the total 
energy enclosed in the sphere at rest mass a. Its change is determined by the surface work, involving pressure p and 
velocity u, and the boundary luminosity L = Aitr 2 q. 

We will now show that Eqs. (§)-(§) can be written in conservative form. Optimal for a numerical implementation is 
a formulation that contains quantities that arc either conserved or local. Since the rest mass as independent variable 
is trivially conserved, we look out for the conservation of volume, energy and momentum. Equation (JsJ) leads to the 
definition of the integral 
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as volume and Eq. (0) suggests the definition of the integral 

r(l + e) 
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as gravitational mass, equivalent to total enclosed energy. For the radial momentum we choose u (1 
show later that it indeed leads to a conservation equation. Borrowing the notation from Romero et al. |lf 
the candidates for conserved quantities and define specific volume, specific energy and specific momentum as 
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In the nonrelativistic limit, we retrieve with a = T = 1 the familiar specific volume D = 1/p, the sum of the specific 
internal, kinetic, and gravitational energy r = e + u 2 /2 — m/r, and the specific radial momentum S = u. With the 
help of Eqs. (||) and @, the time derivative of 1/D leads to the continuity equation (|l^). An energy equation is 
obtained by taking the time derivative of specific total energy 1 + r and substituting Eqs. (|J) and (fjj): 



d 
adt 



T(l + e)+u- 



l__d 

a da 



[47rr 2 a (up + Tq)~^ 



(14) 



This immediately leads to the conservation equation ( |l6| ) for the total energy. A tedious but straightforward calculation 
based on Eqs. (|l4|), (|A3|)-( A10 ), and ( A13 ) finally leads to the momentum equation (|l7|), which completes the set of 
conservation equations 
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(15) 
(16) 
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The time derivative in the last equation usually is small, so that this equation rather acts as a constraint on the lapse 
function a. This equation is derived in the Appendix from the space-space component of the four-divergence of the 
stress-energy tensor. The constraints ( |l8| ) and ( |H^ ) explain themselves in analogy to the Newtonian limit, where the 
first becomes the definition of the rest mass density and the second the Poisson equation for the gravitational potential. 
This set of conservative equations is fundamental for the following discussion. We will first point out their relation 
to the general relativistic jump conditions at a shock front and then derive a consistent incorporation of artificial 
viscosity in Eqs. (|l5|)-(|20|). In a second and third part, we investigate the general relativistic Boltzmann equation 
with two different natural descriptions of the momentum phase space for the radiation particles. Both formulations 
will show interesting relations to the above-mentioned conservation laws. 



III. SHOCK WAVES AND ARTIFICIAL VISCOSITY 



The jump conditions across a shock front reflect the governing conservation laws and can directly be determined 
from the latter: In an isolated shock wave, we assume a stable physical state at the left-hand side (subscript L) and 
at the right-hand side (subscript R). The direction of shock propagation will not play a role in this derivation. A 
conservation law can then be integrated over an infinitesimal range around the shock front, containing rest mass Aa^ 
on the left-hand side and Aan on the right-hand side of the shock at position a s . For example, the continuity equation 
(H) leads to 
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As the shock moves with da s /dt = da^ldt 



dan/dt through rest mass the first jump condition reads 
0. 
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The brackets denote the difference of the enclosed expression evaluated on both sides of the shock. With an analogous 
argument applied to the energy conservation equation, the general relativistic jump condition in spherical symmetry 
with energy transport included reads 
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(22) 



At this point one might guess the correct jump conditions for the momentum conservation from Eq. ( |l7| ) - in spite 
of the complications introduced by the nonvanishing source term. A rigorous derivation along the lines given by May 
and White Ol indeed leads to 



da s 
~dt 
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(23) 



Of course, with q = 0, one immediately recovers the jump conditions for pure hydrodynamics found in Ref. ||. 

Although ideal hydrodynamics allows discontinuities in the physical state across a shock front, there are many 
numerical schemes that cannot handle infinite gradients in the velocity, density, or temperature profiles. A well-known 
solution to this problem is the introduction of an artificial viscosity. The artificial viscosity limits the gradient of the 
profiles and broadens the shock, which then becomes spread over several grid zones. With a conservative formulation, 
the dynamics of these asymptotically incompressible zones becomes determined by conservation of volume, energy, and 
momentum flux across them. This guarantees the conservation of these quantities over the whole shock as required 
by the jump conditions. Older schemes use artificial viscosity in a form proposed by Von Neumann and Richtmyer 
& 

In spherical symmetry, the formulation of artificial viscosity has to be chosen with care in order to avoid systematic 
artificial heating during homologous compression. An excellent solution has been found by Tscharnuter and Winkler 
p2[ for nonrelativistic hydrodynamics. We extend this approach to general relativistic hydrodynamics and define the 
viscous tensor 



Q a0 = Al 2 pu^ 



£a/3 



if < 0, 



= otherwise; 

where P a /3 — u a up + g a /3 is the projection operator onto the three-space orthogonal to the fluid four-velocity u a . The 
artificial viscosity is based on physical viscosity and is chosen according to the standard shear viscosity (Misner et al. 
p0| , exercise (22.6)). It is weighted with a variable coefficient that reflects the local density and compression. The 
length scale Al sets the order of magnitude of the desired shock width. 

First we note that the artificial viscosity tensor is traceless and that in our comoving coordinates its nonvanishing 
components are: 
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(24) 
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This viscous tensor is included in the derivation of the Einstein equ ations in spherical symmetry in the Appendix. 
We can check the behavior of viscous heating in the energy equation (A12): 
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The first expression shows that the viscous heating vanishes in the case of homologous compression (u/r = du/dr). 
This is due to the nonisotropy of the viscous pressure: With homologous compression, the work done on a fluid element 
by radial compression against Q a a would heat the fluid element. But the homologous compression comes together 
with a simultaneous compression in $ and tp direction. By choosing the viscous pressure components to have negative 
sign in these directions, the net heat production vanishes. The second expression together with Eq. ( p4| ) shows that 
viscous heating always has positive sign. In a conservative formulation, the viscosity affects the equation for the total 
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energy evolution ([l6]), the momentum evolution (fi~7|), and the constraint for the lapse function (^). These equations 
then become 



d_ 

m 



dr 


d 


dt ~ 


da 


dS 


d 


~dt ~ 


da 




a 


[( 




r 





1 



47rr 2 p p 



1 



3 (p - Q) 
P 



.da Id, , - , ,i 

^ + ^d-a-- p d-a [aP] -- p d-a [VaQ] 



m 8irr 2 , , . , „, ON 2p Q 

~ + (p 1 + e )(P + Q )-q 2 )-- + - 

r p v ' p p 

1 9 



IV. BOLTZMANN RADIATION TRANSPORT 



A. General relativistic view 



The general relativistic Boltzmann equation in comoving coordinates was derived by Lindquist ||. Radiation that 
is not necessarily in equilibrium with the matter is described by a distribution function /: 

dN = f(x,p) {-p^) drdP. 

The volume element dr is crossed at four-location x by dN radiation particle world lines with four-momenta p in 
the range dP while the observer moves with four- velocity u M . Lindquist || shows that this definition makes the 
distr ibution function Lorentz invariant. We measure the particle four-momentum in a comoving orthonormal frame 
(Al), with the components 



P" 



pcosW, p 







psmu coscp, p 



p sin 6 sin i 



(25) 



The absolute value, p, of the three-momentum can be determined from the scalar particle energy E 



P' 



\/ p 2 + m 2 and the rest mass (the mass of radiation particles is assumed to be zero in this paper) . The direction of 
the particle three momentum is specified by the angle cosine p = cos 9 to the radial direction. In spherical symmetry, 
the distribution function does not depend on the azimuth angle (j). Thus, the particle distribution function depends 
on four arguments 



dN = f{t 1 a,p 1 E)E 2 dEdp, 
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The Boltzmann equation for metric (|l|) then reads (Lindquist 
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(26) 



with A = log(r'/r). The right-hand side of the equation is the collision term that describes changes in the particle 
distribution function due to interactions with matter. It is represented here by an emissivity j and an absorptivity 
X- Since the independent space coordinate is enclosed rest mass, it is convenient to introduce the specific distribution 
function F(t, a, p, E) = f(t, a, p, E) / p. With the help of the relations in the Appendix, the Boltzmann equation (|2T 
can be rewritten in the conservative form (Mezzacappa and Matzner Yamada et al. |3^| ) 
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The first term is the temporal change of the particle distribution function. The second term counts the particles 
that are propagating into or out of an infinitesimal mass shell. The third term corrects for the change in the local 
propagation angle /i when the particle moves to another radius. The curved particle trajectory in general relativity 
is accounted for by the term proportional to the gradient of the gravitational potential 

1 da _ <9$ 
a dr dr 

The term on the second line accounts for angular aberration: When the comoving observers change radius, the 
corresponding relabeling of particle propagation angles also adds a correction to the in-out flow balance for the 
particle distribution function taken at constant /i. The same is true for the frequency shift described in the third line: 
The first two terms in the parentheses account for the Doppler shift caused by comoving observers. The relativistic 
term [iTdQ/dr describes the redshift or blueshift in the particle energy that applies when the particles have a velocity 
component in the radial direction (/i ^ 0) and therefore change their position in the gravitational potential. 

The integration of the Boltzmann equation over momentum space, spanned by the propagation angle and particle 
energy, is supposed to reproduce the conservation laws for particle number and energy as stated in Eqs. ( |15| ) and 
(|l6|). We define H N and G N to represent the zeroth and first /i moments of the distribution function: 

/l />OC 
J FE 2 dEdp, 

"1 />oo 
iN _ I / TPJP2, 



G = J J FE'dE^dp,. 

Integration of Eq. ( p7| ) over momentum space gives the following evolution equation in terms of these moments: 
9H N , d u 2 „ m f j 
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^- Unr 2 apG N ] - a [ J -E 2 dEdn + a ( xFE 2 dEd\x = 0. (28) 
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The aberration and frequency shift terms do not contribute because (l — /i 2 ) vanishes at /i = ±1 and E 3 F is zero 
for E — as well as for E = oo. The nature of Eq. (|2^) is a continuity equation analogous to Eq. jl^) , extended by 
source and absorption terms for the radiation particles. One more integration over the rest mass a from the center 
of the star to its surface finally gives an equation for the total radiation particle number variation as a function of 
integrated emission and absorption. 

Slightly less straightforward is the reproduction of the energy conservation. Defining the energy moments as 



H = J FE A dEdp 

G E = J FE 3 dEpdn = q/p 

P E = [ FE 3 dEp 2 dn, 



the corresponding integration of the Boltzmann equation ( p7j) first leads to the radiation energy equation and the 
radiation momentum equation 

^ + I [4nr 2 apG E ] + ™ (H* - P E ) - + ^ P E 

at oa L v ' \ at r J 

+ T^-G E -a J ] -E 3 dEdp + a J X FE 3 dEd[i = 0, 



dG l 

dt ' da 1 ^ 1 \ r dr J K 1 \ dt ' r J 

+ T^-P E + a X FE 3 dEpdp = 0. 
or J 

Having Eq. (|lj) in mind, we construct, in analogy to T(l + e) + uq/p, the specific radiation energy TH E + uG E and 
investigate its time derivative. After a fair amount of algebra and the use of relations from the Appendix, most of the 
contributions cancel and one is left with the expected result 
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e , U Q E \ + » U^ap (uP E + TG E )] 



aT / -E 3 dEdp + aT / X FE 3 dEdp + au / \FE 3 dEpdp 



(29) 



The time derivative contains the conserved energy; the second term describes the change of energy in the fluid clement 
due to surface work by the radiation pressure, P E , and in-out flow of radiation at luminosity L = Aur 2 pG E . Note 
that (uP E + TG E ), being the first p moment of (TH E + uG E ), indeed is the flux of the conserved quantity at 
the zone boundary. The same feature is realized in the radiation momentum equation (|l7|): d/dt(uH E + TG E ) = 
—d/da\Anr 2 ap(TP E + uG E )] + . . .. It underlines the natural appearance of this formulation. The source terms in 
Eq. (|29| ) describe the energy exchange with matter by emission, absorption, and work due to particle stress. Thus, 
by splitting the internal energy e — e + H E and pressure p = p + P E into contributions from the matter stress-energy 
tensor and the radiation stress-energy tensor, the exact energy conservation equation can finally be written concisely 
as 
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(30) 



B. Order v/c limit 

The 0{v/c) Boltzmann equation was derived by Castor p4] ]. One can also obtain it by directly dropping the higher 
order terms from Eq. (27), i.e., by setting a = T = 1 = const, and replacing u by the nonrelativistic velocity v: 

f + ^[4^] + i|-[(l-^)F] 

, / d In p 3u \ v \ 1 d r „, -, j , . 

The conservation properties of the 0(v/c) Boltzmann equation become apparent in terms of the moments that are 
directly derived from Eq. (|3l|): 

j 



E dEdp + J X FE 4 dEdp = 0, 
X FE 3 dEpdp = 0. 

When we construct the conserved specific radiation energy (TH E + uG E ) in the 0(v/c) limit, we first keep all terms 
that arise from the 0(v/c) Boltzmann equation and obtain 

0=§- t (H E + vG E ) + §-[^p(vP E + G E )] 

3 -E 3 dEdp + J X FE 3 dEdp + v J X FE 3 dEpdp 

Note that the first order terms 
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- r ( HE - pE )-{-^r + v) p 

arising from frequency shift and angular aberration in the radiation energy equation cancel with the zeroth order 
terms in the radiation momentum equation after the latter are multiplied by the fluid velocity. In a numerical 
implementation, the finite differencing of the Boltzmann equation has to guarantee the same cancellation in order to 
achieve 0(v/c) energy conservation. The nonconservative term found in the third line in Eq. (|3^) is of second order 
in v/c. 



V. APPROXIMATE TANGENT RAY DESCRIPTION 



A. General relativistic view 



The description of the particle four-momenta in the comoving frame eases the numerical implementation of the 
collision term. The angle dependent integrals over reaction cross sections of particles with matter are readily evaluated 
if the target is at rest in the comoving frame. On the other hand, the left-hand side of the Boltzmann equation (]27]) has 
to take into account the Lorentz transformation between adjacent comoving observers in the description of the changes 
in the particle distribution function due to propagation. An adequat discretization of the numerous correction terms 
with partial derivatives is, in spherical symmetry, a resolved challenge (Mezzacappa and Bruenn fl35| , Liebendorfer 
|p8|). However, this does not prevent consideration of other interesting options. Here we would like to look into a 
description of the momentum phase space in variables that remain constant along a propagation path in the absence 
of interactions. Particles change their coordinates in the momentum phase space only due to collisions, evaluated 
on the right-hand side of the Boltzmann equation. The left-hand side does not require corrections involving partial 
derivatives with respect to the coordinates of the momentum phase space. The drawback is given by the unavoidable 
need of transformations to relate the cross sections in the fluid rest frame to the distribution function in the chosen 
momentum phase space description. This idea was investigated in polar slicing and the radial gauge by Schinder and 
Bludman |11| . We explore a similar ansatz in orthogonal comoving coordinates in this section. 

A geodesic in a static spherically symmetric space-time can uniquely be described by an impact parameter b and a 
particle energy e at infinity. Let us start with a Schwarzschild metric 

ds 2 = -r|dt| + T s 2 dr 2 + r 2 (dtf + sin 2 Mip 2 ) , (33) 

where we add a subscript S to quantities that could be confused with corresponding quantities in the comoving frame 
(e.g., Ts = y/l — 1m I r). The trajectory of free propagation in the plane of constant ip with energy Es is related to 
the impact parameter b and the energy at infinity e as (Misner et al. fl30| , Eqs. (25.55) and (25.18/19)) 

1 dr\ 2 2m r 2 , 

tm) +1 -- = ¥ (34) 

r s E s = e. (35) 

The relation between the particle propagation in the r- and $ directions is given by the angle 9 the particle trajectory 
makes with the outward radial direction: 



\/l - ri 2 q rdd ,„„. 

tan6> = ^ ^ = — | — . (36) 

r s l dr 

Equations ([Mj) and (|3^) are solved for the relation between the particle angle, fj,s and the impact parameter: 



b=— ^1-/4- (37) 

Our next step is to transform /is and Es into the particle angle, and energy, E, measured by the comoving 
observers. From the metrics (0) and ( |33|) one obtains the relations 

dt s _ aT 



dts_ = 

da r r 



2 (38) 

s 
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for the Lorentz transformation between the Schwarzschild- and comoving-coordinate time. The Lorentz transformation 
also links the particle four-momenta in these two coordinate bases: 



Ps= ry^s/i^-V 1 -^ ) E s 



V> t =(a- i ,^ii,^y/x=tf,0)E. (39) 



From the application of transformation (38) to the four- momenta in Eqs. ([39|), we extract the transformation of the 
particle propagation angle and energy: 

1-/4 _ 1-m 2 



r l (r + u/ i) 2 

T S E S = (T + un) E. 
Finally, Eqs. (p7|) and (pq) become, in terms of comoving frame variables, 



> = r£HZ (40) 

1 + UfJ. 

e = (r + it/x) E. (41) 

These useful relations link the local particle angles and energies measured by comoving observers at different radii and 
with different velocities to the particle impact parameter and energy at infinity. Additionally, these relations provide 
insight into the nature of the energy conservation equation (E9^. The conserved quantity under the time derivative 



TH E + uG E = J (r + UfJ.) FE 3 dEd^ = J eFE 2 dEd^ 



is the total radiation energy at infinity expressed in terms of the comoving frame radiation energy and momentum. 
This also makes clear the following important point: The locally observed radiation quantities must be transformed 
to a common observation point (e.g., infinity) before they can be integrated to define a conserved quantity. 

Let us recall that the particle trajectories along constant impact parameter, 6, and constant energy at infinity, e, 
were derived in static vacuum Schwarzschild space-time surrounding a gravitational mass m. Nevertheless, wc can 
span the momentum phase space in the dynamical Boltzmann equation in comoving coordinates, Eq. (|27j), with 
the coordinates Q (b,e) instead of the comoving frame four-momenta (/i,E). If the static trajectories are a good 
approximation to the trajectories on a dynamical background, we expect that the momentum state of free particle 
propagation, measured in (b,e), barely changes between adjacent comoving observers. This would imply that the 
terms involving partial derivatives with respect to b and e become small compared to the time derivative of the 
distribution function and the in-out flow term on the left-hand side of the exact Boltzmann equation. 

We replace the partial derivatives in the Boltzmann equation by directional derivatives along the static particle 
trajectory: We add an overbar to the distribution function f(t,a,b,e) — f(t,a,[i,E) in order to indicate that the 
partial derivatives are taken at constant impact parameter b and energy at infinity e and make the following ansatz: 



The coefficients Cj are calculated by comparing the Boltzmann equation (|26|) to Eq. (|42|). The relations (J40j) and 
( [4l| ) are used in a straightforward but lengthy replacement of the partial derivatives in Eq. (^6|) by derivatives with 
respect to impact parameter and energy at infinity. The Boltzmann equation written in terms of these directional 
derivatives reads 



1 Strictly speaking, the pair (6, e) only specifies a trajectory. The propagation direction on this trajectory can, for example, be 
addressed with the convention that a negative impact parameter denotes a propagation with decreasing radius and a positive 
impact parameter a propagation with increasing radius on the same trajectory. 
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10/ r 0/ 4irr / / p\ . a . 

a of r aa 1 + up \ \ p J ' 

Its left-hand side (LHS) involves the first two ususal propagation terms that relate the change in the distribution 
function to particle in-out flow at the boundaries of a mass element. The third term arises owing to the drift of 
the particle location in the phase space from constant (6, e) in front of a non-vacuum background. This general 
relativistic term is of order (Gpr 2 /c 2 ) and vanishes in the vacuum. It is interesting to see in the static limit that 
this term mainly arises from the background presence of matter and not from its dynamical motion. We might then 
neglect this 0(Gpr 2 / c 2 ) term in the exact equation ([43j ) and find to the very simple approximative general relativistic 
Boltzmann equation in spherical symmetry: 

Id] Yd] 

-777 +M-Tp = J +XJ- (44) 
a at r da 

This approximation is excellent in the proximity of compact objects; it includes full gravitational redshift. It might 
however fail within regions of extremely high energy density. Moreover, we show in the next section that its 0(v/c) 



expansion is identical with the full 0(v/c) Boltzmann equation (31). Equation ( |44| ) therefore also incorporates an 
angular aberration and Doppler shift between the comoving coordinate frames. Its intuitive form simply reflects 
Lindquist's general starting point B (2.12/22), where the Boltzmann equation is written in terms of directional 
derivatives of the distribution function along the phase flow. 



B. Order v/c limit 



First, we multiply the left hand side of the 0(v/c) Boltzmann equation ( |3l| ) by the density, p, express it in terms 
of the neutrino distribution function / = pF, and take the prefactors out of the derivatives to obtain the formulation 
presented by Castor l34j . With the continuity equation 



d lnp 

dt 



2v 
r 



we then find 



LHS 



0/ 
dt 



0/ 
Of 1 



E 



0/ 
dE' 



(45) 



(46) 



On the other hand, we expand the approximate general relativistic Boltzmann equation ( |44| ) to order v/c and 
compare to Eq. ([46|). Before the expansion, however, we have to transform the directional derivatives back to partial 
derivatives with respect to comoving frame energy, E, and angle, p. In terms of an unspecified parameter A, Eqs. 
( |40| ) and (^l|) lead to the derivatives of the comoving frame propagation angle and particle energy, taken at constant 
impact parameter and energy at infinity: 



Op 
0A 

dE 
0A 



(1 



1 



b,E 



"AO 
E 

+ M 



vp 



(1 



1 dr 
r dX 

2\ v ® r 
' rdX 



dv 



\- vp dX 
, dv 



dX 



(47) 



The evaluation of the partial derivatives on the left hand side of Eq. ( |44| ) is straightforward with the relations (|4 
We obtain 



= dj_ 
dt 



1 



d]_dp_ 

dp dt 

r' da 



d]_dE_ 

dE dt 



1 



p df 
- — — 4 

r' da 

(l + vp) 



^d£dp 
r' dp da 
1 p 



r 1 + v I p ' 



p_d]_dE_ 

r' dE da 
df 
dp 



r 1 + v I p r' 



E 



0/ 
dE' 



(48) 
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Following Bruenn p3], we have already dropped the terms involving the matter acceleration dv/dt in the last step. If 
we further neglect the order v/c and higher order terms, Eq. ( E ) reduces exactly to Eq. ([46|). We therefore conclude 
that the approximate general relativistic Boltzmann equation (j44j) exactly includes and, moreover, extends the 0(v/c) 
Boltzmann equation in the presence of strong gravitational fields. 



VI. CONCLUSION 



We provide the exact Einstein equations for radiation hydrodynamics in spherical symmetry, formulated in the 
comoving orthogonal coordinates introduced by Misner and Sharp |§,^6|. We show that these equations can be 
written in a concise and strictly conservative form. The general relativistic jump conditions at shock fronts are 
derived under inclusion of the radiation energy flux and radiation momentum flux. 

The conservative formulation is ideally suited for the application of numerous numerical schemes specifically designed 
for the solution of partial differential equations in this form: for example, adaptive grid techniques (Winkler et al. 
p7|), Riemann solvers (Godunov |^8| , Davis |27|| ), and high-resolution shock capturing schemes (Romero et al. 
Marti and Miiller p9| ). Some of these schemes require artificial viscosity for a numerical representation of shock 
waves. We have generalized the tensor artificial viscosity ansatz of Tscharnuter and Winkler |32j to general relativity. 
As in the nonrelativistic case, viscous heating is shown to always have positive sign and to vanish during homologous 
collapse. 

In the second part of this paper, we focus on the general relativistic Boltzmann equation in our chosen coordinates. 
We confirm that, as expected, the appropriate moments of the particle distribution function, evolved according to the 
Boltzmann transport equation, obey the conservation laws found in the first part of the paper. In the 0(v/c) limit, we 
identify leading terms that require careful discretization in a finite difference representation of the Boltzmann equation. 
Angular aberration and frequency shift corrections arise in the propagation part of the Boltzmann equation since the 
time and space derivative of the particle distribution function is taken at constant comoving frame four-momenta in 
the particle momentum phase space. 

In the last part, we try to recover the simple nature of the left-hand side of the transport equation, understanding it 
as a directional derivative along the phase flow of the particles. The approximation of the phase flow by geodesies in 
vacuum Schwarzschild space-time leads first to the exact general relativistic Boltzmann equation with the momentum 
phase space now described in terms of particle impact parameter and energy at infinity. The correction terms to the 
partial derivatives with respect to time and space along vaccum geodesies are of order (Gpr 2 /c 2 ) in this representation. 
They are only due to the background matter distribution neglected in the determination of the phase flow. We 
show that neglecting these terms leads to a very intuitive Boltzmann equation that is exact in the 0(v/c) limit and 
additionally includes gravitational redshift. Although the numerical implementation of this equation has to solve other 
problems (for example, the radius dependent range of valid impact parameters b — [r, oo]), the exploration of this 
idea reveals the physical interpretation of the conservation laws found in the first and second parts: The conserved 
quantities are those measured by an observer at infinity, but expressed in terms of local quantities measured by 
comoving observers. 

After Wilson's pioneering work, years of code development by several groups (Mezzacappa and Bruenn ||35"| , 
Yamada et al. 33 1, Burrows et al. Liebendorfer pii) ], Rampp 42 1, Messer ^3|) have led to Boltzmann solvers for 
neutrino transport in spherical symmetry. Dynamical simulations of core collapse supernovae have been presented in 
the Newtonian limit by Mezzacappa et al. Q] and Rampp and Janka [E5[. As general relativistic effects significantly 
influence the postbounce evolution of a collapsed star (Bruenn et al. J46|X these codes have ultimately to be extended 
to include full general relativistic radiation hydrodynamics. Such simulations were performed by Liebendorfer et al. 
p8| , p9| based on the equations discussed in the first and second parts, under omission of neutrino back reaction. 
Awareness of the conservation laws is essential in the simulation of supernovae because the observed explosion energy 
is two orders of magnitude smaller than, for example, the core binding energy or the radiation (neutrino) energy. It 
has to be excluded that an energy drift, owing to numerical inaccuracy in energy conservation, seriously affects the 
equalized balance and causes or prevents an explosion. 
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APPENDIX A: DERIVATION OF THE EINSTEIN EQUATIONS IN SPHERICAL SYMMETRY 

In this appendix, we derive the Einstein equations in comoving coordinates and spherical symmetry. We follow 
closely the guidelines of exercise (14.16) in Misner et al. |Q. In spherical symmetry, we are allowed to make the 
following ansatz for the metric: 



2 i -2 

- 1 sin 



The coordinates consist of an independent time coordinate t and a space coordinate a. The areal radius r(a, t) is chosen 
such that the area of the two-sphere, dd 2 + sin 2 (i?) dip 2 , is Airr 2 . We define the following orthonormal noncoordinatc 
basis 

w * = e^dt 
u a = e^^da 
w* = r{a,t)dd 
Ld v = r{a, t) sm(i))dip 

and calculate its exterior derivatives. From = c?w M + uj^ A u> v (14. 31a, b in Misner et al. 
nonvanishing connections 

— <E» ,a 



(Al) 

J) one determines the 



UJta = 


-$'e~V - 


u*e = 


r 


&tip — 


— e 
r 


u a § = 


r 


l^aip — 


r ' -A m 

e 

r 


k>&<p — 


1 cos(tf) 

— ■ CtJ 

r sin(i9) 



(A2) 

An overdot denotes the derivative with respect to coordinate time t and a prime denotes the derivative with respect 
to the spatial coordinate a. The exterior derivatives of these connections lead to the curvature tensor, and finally the 
Einstein tensor: 



(F + 2F) 



2H 



G a 

G\ = -{F + 2E) 
G% = - (E + E + F) 
G v n = -(E + F + E) 



where 



E 
E 
H 
F 
F 



(_$" 

l r 



e -2A 



r 

- I r 

r 



$'A' - $ 
(r - r$) e" 2 * - r'$'e~ 
r$' - r'A] e -(*+ A ) 



A - $A + A 1 



4(l + f 2 e- 2 *-r' 2 e- 2A ) 



' a r\ „-2A 
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In the following, we try to write the Einstein equations G = 8ttT as concise as possible in terms of quantities that 
are defined in the restframe of a fluid. The nonvanishing components of the stress-energy tensor of a fluid in its 
comoving frame are given in Eq. (|J). We include an addition, Q, to the diagonal component, representing artificial 
viscosity as defined in Sec. Ill, Eq. (^4|): 



p 



T 







P 



First, we eliminate the exponentials in the metric by substituting the lapse function a for e* and the function r'/T 
for e A in order to retrieve notation (|l]). Then, we have to make sure that the yet unspecified space coordinate a is 
attached to matter, i.e., that it is a Lagrangian (comoving) coordinate. The rest mass between coordinate ag and a\ 
at a fixed coordinate time t is given from the metric by 



A(a ,ai) 



2tt />tt />ai 
Ja 



r' f 
p—da rdd r sin(i9)dip = I 

* J an 



4irr 2 r' p 



da. 



We tie the spatial coordinate to rest mass by requiring that A(ao, ai) = da for arbitrary boundaries ao, a%. This 
is equivalent to adopting the relation 



r = 



Airr 2 p' 



(A3) 



Further, we define a "velocity" u — r/a that describes the change of areal radius with proper time of the comoving 
observer. 

From the nondiagonal component of the Einstein equations we derive in three steps: 




~v! 1 (H 




r' a \r' 


-f)l 






r ' 




' Aitrq 




~ + T • 





47rg, 



(A4) 
(A5) 



Equations ( |A3| ) and (A4) lead in three further steps to the evolution equation for the rest mass density: 

d fl\ d fATtr 2 r r 
adt \p 




(A6) 



Equation (A4) can also be combined with the time-time component of the Einstein equation to give in three more 
steps: 



1 



G 



X -r 2 r'G tt = r -(l + u 2 -T 2 )+rr' 



V 


--( 






a r \ 


K r> - r) 


TV 


uv! 


Airruq 


— 




r 



-(1 

2 V 



-T 2 ) = Anr 2 r'p (1 + e) + 



Aitr 2 r'uq 



(A7) 



It is now convenient to define a "gravitational mass" m = r/2 (l + u 2 — T 2 ) that leads, with the help of relation ( [A3] ) 
and (A7), to the concise expressions 
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r=^l + u 2 -^ (A8) 

m' = r(l + e) + ^. (A9) 
P 

The physical interpretation of m(a,t) as the total energy of the fluid inside coordinate a becomes clear if one looks at 
an evolution equation for the mass m. In order to der ive such an equation, we take the time derivative of Eq. QAS| ). 
We already have the time derivative of T in Eq. ( |A5| ) and the time derivative of r in the definition of the velocity, 
but we still need to calculate the time derivative of the velocity before being able to isolate the time derivative of m. 
An evolution equation for the velocity can be deduced from the space-space component of the Einstein equations: 

1 ™„ „ ^ 1 m 1 fii ^ 2 $' 



G aa = -F-2E=-- 1--T 2 -) =Att( P + Q), 

2 r r r \a r 1 



-=T 2 %-™-4Trr(p + Q). (A10) 

Therefore, the evolution of the gravitational mass is 

!^l = -4 7T r 2 (u(p + Q)+r q ). (All) 

The change of the total energy m within a sphere is given by surface work against the pressure and the in-out flow of 
energy by transport processes. 



Analogous to th e der ivation of the density evolution equation ( A6) it is possible to take the time derivative of Eq. 
(|A|) and use Eq. ( |All| ) together with one of the angular components of the Einstein equations to derive an evolution 
equation for the specific energy. However, we choose to derive these equations more simply from the vanish ing four- 



divergence of the stress-energy tensor. First, we list out the nonvanishing connection coefficients from Eqs. (A2) and 
the definition w M „ = T ^ ua u a that are used in the following computation of the four-divergence: 
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Note that we have to respect T a ^ in our noncoordinate basis (Al). With this in mind, one easily obtains an 



evolution equation for the specific energy from the time component of the four-divergence: 

= -T% 

P ' 

de i ^ d { 1\ 1 fA _ 2 2\' 'AuQ 



-adt' {p + Q) ^dt\- p ) + ^^ r ^-Vf (A12) 

The last equation - usually used to determine $ and the lapse function a = e* - is then derived from the space 
component of the four-divergence: 

r' 1 

ri rpCLU 

Yp ' v 

= + ( 1 + e + 2±R)# + A. f_i_!) + *LQ. (A13) 

P \ p J ceot \ATrr z p p ) r p 

This set of equations (without the artificial viscosity) was, with some approximations, derived by Misner and Sharp 
~36l. 
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